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THE APPROXIMATE SOLUTION OF SINGULAR INTEGRAL EQUATIONS 
ARISING IN ENGINEERING PRACTICE 


Francis B. HILDEBRAND 


Presented Feb. 7, 1941, by E. C. Kemble and J. H. Van Vleck. 


1. Introduction. In a recent paper! Crout has 
given a method of solution of integral equations 
by polynomial approximation. The purpose of 
the present paper is to extend this method to the 
solution of certain integral equations which arise 
in a natural way in the formulation of problems 
of static field theory, air-foil theory, elasticity, 
heat and fluid flow, diffusion, etc., wherein the 
unknown function, as well as the kernel or 
Green’s function, may have infinite singularities 
of known order. 

Several specific problems of this type are solved 
numerically, the problems being chosen, for the 
most part, so that the exact solutions are known, 
and the results are shown to be well within 
engineering limits of accuracy. Use is made of a 
method of least squares which is given in a paper? 
now in publication, and somewhat extended in 
the present paper. 

2. Outline of the General Method. We consider 
here the solution of integral equations of the 
first kind, of the form 


gx) 


where 9 and G are given, and suppose that o(2) 
is known, from physical considerations or other- 
wise, to have singularities of known order at 
certain points in the interval (a, 6). Constructing 
s(x), a linear combination of r suitable functions 
s(x), possessing the correct singularities at these 
points, 


s(2) = 


where the y’s are constants to be determined, 
we set 


(2.2) o(x) = a(x) = s(x) + y(2), 


where y(x) is a polynomial of degree 2n, specified 
by its ordinates at 2n + 1 points equally spaced 
over the interval (a, b). 

If the spacing is h and the corresponding 
values of y(x) are denoted by y—n, Y-nit, 5 
Yor. 5 Yn—1y Yn, and if we set x = (a+b)/2 
Ah, the polynomial can be expressed conveniently 
in the Lagrangean form 
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= 


=—n 


(2.3) 
where K,() is a polynomial of degree 2n in i, 


2n 
(2.4) Kia) = Baik’. 
s=0 
Tabulated values of the A’s for —n < % < n by 
tenths, n = 1, 2, 3, as well as tables of the co- 
efficients B,;, are given in Crout’s paper.' 
Now G(x) is specified by 2n + r+ 1 indepen- 
dent parameters. Substituting ¢(x) for ofa) in 
the right hand side of (2.1) and assuming that 
approximates satisfactorily, it follows 
that 


r b 
(2.5) ox) = + 
j =] n b 

i=-—n @ 
The parameters can now be determined by any 
set of conditions leading to a set of 2n + r+ 1 
independent linear equations. Suppose that we 
require the integral of the square of the difference 


between the two sides of (2.5) to be a minimum.’ 
If we write (2.5) in the form 


j= 


t= —% 


the y’s and y’s are to be determined so that 


(2.6) LiS (x) +> fz) | dx 
= min. 


Equating the partial derivatives to zero, we obtain 
the system of linear equations 


i=—n 

(2.7) S x)8,(x)dx, p= 2, 
j=1 


S Y(x)de, . « «98: 


r 
~/ 
r b b 
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In general, direct integration would be im- 
practical. However, if we employ approximate 
integration with weighting coefficients D,, this 
system of equations takes the form 


i=—yn 


k=1 


(2.8) 
(ax) ¥ ots) + (are) (xx) 


=] 


k=1 


It is now easy to verify the, fact that the aug- 
mented matrix of equations (2.8) can be obtained by 
multiplying the augmented matrix of the system 


(2.9) 9(ax) = + (ae), 


i=—n 
.»m, 


by the matrix formed by multiplying the kth column 
of the transposed coefficient matrix by the integra- 
tion coefficient D;,. For, if we denote S;(a;,) by 
S* and Y;(a,%) by Y,*, this matrix product takes 
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the numbers D;’ = uD;, where uv is chosen con- 
veniently. This gives a practical method of ob- 
taining 2n + r+ 1 linear equations to determine 
the parameters specifying the approximate solu- 
tion (x). In each case a rough check on the error 
ean be had by substituting the calculated values 
of the parameters in equations (2.9). 

The examples which follow will illustrate the 
fact that the coefficients of (2.9) can be obtained 
conveniently and systematically by the use of 
matrix multiplication, wherein certain matrices 
depend only upon the particular kernel involved, 
while others which depend only upon certain 
general characteristics of the problem can be 
computed once and used repeatedly in a variety 
of problems of a similar type. 

3. Applications. To illustrate the procedure 
we first set the following problem: 


Problem 1. A long, uncharged, perfectly con- 
ducting strip of breadth 2a is placed in a two- 
dimensional electrostatic field (a field which does 
not vary in the direction of the length of the 
strip). Find the distribution of the induced charge. 
We consider a unit length of the strip and field 
and denote by o(x) the required charge density, z 
being measured from one edge of the strip along a 
perpendicular to the edge. Then, choosing the 


the form line x = a, at the center of the strip, to be the 

DS? DS? 
DeY,2....Dm¥n" || 


and, multiplying according to the row-times- 


locus of zero potential, the potential at a distance 


column convention of matrix multiplication, the 
product matrix, 


x from the edge 2 = 0, due to a charge distribu- 
tion o(2x) is given by 


JS log 


™m 
where >, has been written for >.>, is seen to be “ld 
k=1 


the augmented matrix of equations (2.8). 
It is also clear that in place of the coefficients 
D,, it is permissible to use as weighting coefficients 


If we denote the sihiadhdae at a due to the exter- 
nal field by —¢(a), again taking the line x = a to 
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be the reference locus, so that ¢(a) = 0, and if we 
assume that the presence of the strip has a 
negligible effect upon the distribution of charges 
producing the external field, then the condition 
that the strip be an equipotential surface is given 
by the equation 


(31) g(x) = og 


Since o(2) is known, this is an integral equation 
of the first kind for determining o(z). 

It is well known that at a non-reentrant corner 
of angle opening 8 the charge distribution on a 
conductor has a singularity of order (1 — =,9). 
In particular, for a strip of finite breadth we have 
§ = 2x at the edges so that the singularities at 
the edges are of order 14; that is, o(x) ~ Aa near 
z=0, and o(x) ~ B(Qa—z2) near x = 2a, 
where A and B are constants. 

It is convenient to solve Problem 1 for two 
special types of fields and to obtain the general 
solution by superposition. 

a. Suppose that the external field 1s antisym- 
metrical with respect to a plane perpendicular to 
the strip at its center line. Then g(a) = 


—o(2a— 2). Here we set s(x) = 


(: _ ‘) where y is an arbitrary parameter, 
a 


and approximate by the combination of 
and a fourth degree polynomial specified by its 
values Yo, 1, Yo, —Y1,and —Yy at the points z = 0, 
a/2, a, 3a/2, and 2a, respectively. Since o(2) is 
continuous at z = a, it follows that yo = 0, and the 
polynomial can be expressed in the Lagrangean 
form 


2 
y(x) = = (A)Yis 


i=-—2 


where A = 2(a/a—1) and K,(A) = K_,(A) — 
Assuming that is approximated 
satisfactorily by o(a) = s(x) + y(x), it follows 
from (3.1) that 


We next set 7 
(3.3) gp(x) = -2f log ds 
GA) gale) = 


and require that 
(3.5) = + = 1, ..., 5, 


corresponding to equations (2.9) of the preceding 
section, where the points 2% are given by x = 0, 
a/8, a/4, a/2, and 3a/4, respectively. Then also, 
because of the antisymmetry, (3.2) is true for 
the corresponding points in the interval (a, 2a) 
and at x = a. 

From (3.3) we have 


If we set 


C. 


and notice that h = a/2, the coefficient of y; in 
the kth equation of (3.6) enema 


so that the matrix of coefficients of the y’s in 
(3.6) is given by the 


where the row and column indices are / and s in 
the first matrix and s and 7 in the second, re- 
spectively. 

The coefficients c,; are found by first placing 
A = 2(&/a — 1) = &/h — 2 in (2.4), so that 


4 
= 2, BalS/h 


Zk 


Csi 


log 


b 


the values of the B’s being taken from Crout’s 
tables.! Then, placing K,(A) = K_(A) — Kid), 
the coefficient ¢,; is read directly as the coefficient 
of (&/h)* in the expression for K,(A). 


(3.2) —2 S| + y(2)} log dé. Carrying out the integrations directly, the 
é—a required matrix is given by the product** 
1.38629 4.77259 12.72690 35.51404 105.11160 6 0 1.11111 1.7777 
91871 4.22212 11.65644 32.71051 96.75405 1 —ll 16 1.23431 2.30078 
.63275 3.61684 10.42869 29.58000 87.57357)| - - 6 .99886 2.70495 
26162 2.34721 7.47519 22.00422 65.80978), © 1 2 46753 2.71285 
06317 1.10524 3.92619 12.21317 37.47033 0 O .14486 1.66806 
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Hence, writing 9p; for gp(2;,), equations (3.6) 
have been obtained in the form 
gp, = 1.11111 ayo + 1.77778 
Qp2 = 1.23431 aye + 2.30078 ay, 


(3.7) ops = .99886 aye + 2.70495 ay, 
gps = .46753 aye + 2.71285 ay, 
= .14486 1.66806 ay}. 


Next, computing 


{ (7) -(2-2) 


directly, we obtain 
Qs1 = 7.84207 ya, ose = 6.56236 ya, 
(3.8) Oss = 5.50583 ya, 9s4 = 3.58146 
= 1.76929 74, 


d 
s—a 


where = 9s(2x). 


Combining (3.7) and (3.8), we have determined 
equations (3.5) in the form 

91 = 1.11111 aye + 1.77778 + 7.84207 ay, 

go = 1.23431 ayo+ 2.30078 ay: + 6.56236 ay, 


(3.9) = .99886 ayo+ 2.70495 ay; + 5.50583 ay, 
= .46753 2.71285 ay, + 3.58146 ay, 
$5= .14486 ayo+ 1.66806 ay; + 1.76929 ay, 


where o = 9(.;). 


In order to apply the Jeast squares method dis- 
cussed above, it is necessary to choose weighting 
coefficients D,’ associated with the points 2. 
Because of the unequal spacing of these points 
the ordinary Cotes integration coefficients can 
not be used. Instead we here require that D; be 
proportional to 2,,, — 2,_,, where we define 


ro = 0, a6 = a,5 and set, for convenience, D;’ = 
2 / 
— 2,1) so that D,’=4, D,’ =, 


D;’ = 3% and D,’ = D;’ = 1 

Now in accordance with the rule given in 
connection with (2.9), the augmented matrix of 
the required set of three equations in ye, yi, and 
y, is obtained by multiplying the augmented 
matrix of (3.9) by the matrix formed by multi- 
plying the kth column of the transposed coefficient 
matrix by D,’, viz., by the matrix 


27778 ~=©.61716 46753 ~—.14486 
44444 1.15039 2.02871 2.71285 1.66806)|- 
1.96052 3.28118 4.12937 3.58146 1.76929 


In this way we obtain the equations 
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(3.10) 
2.05827 ay2 + 5.45016 ay; + 12.28378 ay = 
2777891 + .6171692 + .7491593 + 
4675394 + .1448695, 
5.45016 ays + 19.06646 ay, + 34.87161 ay = 
444449, + 1.150399. + 2.028713 + 
2.71285¢4 + 1.66806¢5, 


12.28378 ay2 + 34.87161 ay; + 75.59963 ay = 
1.960529, + 3.2811892 + 4.1293793 + 
3.9814694 + 1.76929¢5. 


Solving for yo, y: and y,° we have the general 
solution to problem la in the form 


(3.11) 
ay. = — 1.586279, + 1.23571¢2 + 1.86635¢3 — 
8042194 — 1.73238¢6, 


— .269159, — .0322992 + .1767693 + 
.29758¢4 + .16028¢;, 


Cy = 4078391 .142499 — .3301 793 + 0407894 
+ 2309695. 


As an application, we investigate the distribu- 
tion of charge induced on a long, uncharged, 
perfectly conducting strip which is placed in a 
uniform electrostatic field of intensity FE, per- 
pendicular to the center line of the strip and 
parallel to its surface. In this case 9(a) = (a—2a)E, 
and, substituting 9, = aE, = .875ak, 93 = 
dak, and gs = .25aE in (3.11), 
we obtain ye = .05954 FE, y; = .02403 E, y= 
.11366 E. A table of values of the approximate 
solution c(x) is easily computed with the aid of 
the Lagrangean interpolation coefficients tabu- 
lated in Crout’s paper.! 

The exact solution is known to be 


a—x 
= E. 
x(2a—2x) 


In the following table we compare the exact 
values o(x) with the computed approximation 
¢(x) over a half-width of the strip. 


g(x)/E o(x)/E (o—a)/o 
.05 - 48204 -48422 .0045 
.32792 32861 .0021 
21265 21221 — .0021 
3 15655 . 15600 — .0035 
4 .11970 .11937 — .0028 
.09196 .09189 — .0008 
.06952 .06946 — .0009 
a -04985 .05005 .0040 
8 .03229 .03249 .0062 
9 .01588 .01600 .0075 

1.0 -00000 .00000 -0000 _ 
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Near 2 = 0, we find from the exact solution 


Comparison of the coefficient with the computed 
value of y indicates that the “limiting”’ error at 
the edge of the strip is less than 1%. 

b. Suppose that the external field 1s symmetrical 
with respect to a plane perpendicular to the strip at 
its center line. It follows that 9(x) = 9(2a—z2). 
In this case we set s(x) = x¥{(x/a)?+ (2 —- 
z/a)*} and approximate o(a) by the combination 
of s(x) and a fourth degree polynomial y(2) 
specified by its values y2, yi, Yo, 1, and Yo, at the 
points « = 0, a/2, a, 3a/2, and 2a, respectively. 
The polynomial can now be written in the form 


2= 
y(x) = 2, 


where A = 2(a7/a— 1) and = A,(A) + 


K_(X). Proceeding as in part a, the equations 


= 2S is(5) + y(5)} log dé, 
"keh, . oh 


are obtained in the form 


(3.12) 
91 = .04648 aye — .54939 ay, — 
2.76769 ayo — 1.58216 ay, 


= .67112 aye + .42990 ay: — 
2.93845 ayo — 1.01711 ay, 
= 48084 aye + 1.09197 ay, — 
2.83832 ayo — .69117 ay, 
% = .14615 ayo + 1.14969 ay; — 
1.81909 ayo — .28149 ay, 
% = .02020 aye + .40336 ay; — 


54990 ayy — .06747 ay. 


It is obvious from physical considerations that 
the induced charge distribution is indeterminate 
unless the total charge on a unit length of the 
strip is specified. Hence, if or is the total charge 
per unit length, we have the additional condition 
that 


or = + 
or 


(3.13) op = 31111 aye + 1.42222 ay, + 
.26667 ayo + 5.65685 ay. 


Using (3.13) to arbitrarily eliminate y from 
(3.12), applying the method of least squares of 
Problem la to the resultant equations, and 


solving the three equations so obtained for yo, 
yi, and yo, we obtain the solution to Problem 1b 
in the form 


(3.14) 


ay, = —2.713769; + 3.08414q2 + 2.12828¢3 — 
3.58497¢5 — 2.00920¢5 — .14679¢r, 


ay; = —1.02359¢, + .4592992 + .67631o3 — 
.19349¢4 — .2140495 — .13325¢7, 
ayo = —.944629, + .6578892 + .5215393 — 


8584094 — .4918895 — .13077s7, 


ay = 451119; — .31610g2 — .31167¢3 + 
28627¢4 + .18750¢5 + .22452c7, 


the last equation being obtained from (3.13). 

As an application, suppose that a charge of or 
per unit length is placed on a long, perfectly 
conducting strip, no external field being present. 
We require the steady state charge distribution. 
In this case 9(x) = 0, and we obtain from (3.14) 
the values = —.14679s7, ay; = —.133257, 
ayo = — .13077¢7, and = .22452¢r. 

The exact solution is known to be 


oT 
In the following table the exact values o(a) and 


the computed approximate values ¢(x) are com- 
pared over a half-width of the strip. 


= 


\aé(x)/or \ac(x)/or v/a \aé(xr)/or \ac(xr)/or 


.05 {1.02033 |1.01941 6 .34736 | .34730 
1 .73037 | .73025 of .33371 | .33368 
.93023 | .53052 8 .32487 | .32487 
44551 | .44572 9 .31988 | .31991 
4 .39782 | .39789 || 1.0 .31826 | .31831 
.36758 | .36755 


It is seen that all errors involved are less than 
one tenth of one percent of the true values. 

Referring to the exact solution, it is seen that, 
near 2 = 0, o(x) = .22508 (cr/a)(x/a)-, so that 
the “limiting error’ at the edge is less than one- 
fourth of one percent. 

ce. Since any continuous function defined over 
(0, 2a) can be resolved into the sum of two 
continuous functions, one of which is symmetri- 
cal, the other antisymmetrical, about x = a, it 
follows by superposition that we can use (3.11) 
and (3.14) to obtain the approximate distribution 
of charge induced by an arbitrary field of the 
type described in the statement of the problem. 
Then, by the use of further superposition and 
approximate integration, the potential or field 
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intensity at an arbitrary point in space can be 
readily determined. In this sense we have ob- 
tained an approximate solution to the two- 
dimensional Dirichlet problem for an infinite 
strip. 

The integral equation (3.1) can be obtained 
formally by integrating both sides of the equation 


(3.15) 
x 


E(2) = S 


, 
where Ek = - , solutions to which have been given 


under certain restrictions by Betz,’ Fuchs,® 
Schréder,? and others. We have deliberately 
chosen this equation because of this fact, per- 
mitting a verification of the approximate results 
obtained. However, it is seen that the methods 
given here are applicable to very large classes 
of linear integral equations, whether of first or 
second type, and also to many integro-differential 
equations. 

Much of the data tabulated in the solution of 
Problem 1 can .be used in connection with the 
solution of other Dirichlet or Neumann type 
problems associated with various plane figures 
bounded by straight lines. For example, in order 
to determine the charge distribution on a parallel 
plate condenser constructed of long plates of 
finite breadth 2a and separation / it is necessary 
to solve the integral equation!” 


2a = 7 
(3.16) = —2 o(§) log d + 
0 cma 
2a 2 
=) log 
f 
where o(2) = and = o7, the total 


charge per unit length on each plate. An ap- 
proximate solution where h = a/2 has been ob- 
tained by the writer using the methods of this 
paper.’ The approximate charge distribution 
curve for a half plate is given in Fig. 1. Using 
this solution it is not difficult to obtain the 
capacity of the condenser or the potential or field 
intensity at an arbitrary point in space. 

As a final application we set the following 
problem from the theory of elasticity: 

Problem 2. A long bar of equilateral triangular 
cross section 1s twisted by couples applied at the 
ends. Determine the stress function and investigate 
the shear stress distribution in a cross section of the 
bar. It is well known" that the stress function ® 
satisfies the differential equation V2 = —2G0, 
where VY? is the two-dimensional Laplace operator, 
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Fig. 1. 


and vanishes at the boundary of the cross- 
section. Gand @ are constants, G being the mod- 
ulus of rigidity and © the angle of twist per unit 
length of the bar. If we represent the stress 
function in the form 


(3.16) 


where r is measured from the axis of the bar in a 
plane perpendicular to the axis, it is necessary 
to determine a function 9g, satisfying Laplace’s 
equation, such that, along the boundary of a 
cross section, 


GO 


(3.17) 


It is seen that the determination of 9 is analo- 
gous to the determination of the electrostatic 
potential inside the boundary formed by three 
infinite conducting strips coincident with the 
faces of the bar, if the system is placed in an 
electrostatic field such that the potential at a 
point inside or on the boundary, due to this 
field, is given by —GOr*/2. It will be convenient 
to preserve this analogy in our solution to the 
problem. Thus, as in Problem 1, we determine a 
distribution of (fictitious) induced charge on the 
strips necessary to maintain the boundary at 
constant potential. Then, by superposition and 
approximate integration, we are able to obtain 


> 

| 

i 
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the value of 9 and hence, by (3.16), of the stress 
function, at any point inside the boundary. 

It is clear that, in obtaining the potential 
function for the interior of the triangular bound- 
ary, it is necessary to specify the total fictitious 
charge per unit length on the bounding strips, 
the specification being, however, perfectly arbi- 
trary since the potential inside the boundary is 
unchanged if to the original distribution is added 
any multiple of the charge distribution which 
would maintain the conducting surface at zero 
potential in the absence of an external field. For 
convenience we shall require that the total charge 
be zero. 

From symmetry it is obvious that the charge 
distributions on the three faces will be identical 
and, on each face, will be symmetrical about the 
center line. Thus, if we denote the perpendicular 
distance from an edge on any face by the variable 
x, and consider a unit length of the system, we 
are to solve the integral equation 


— (2) -2 log ——— §) d= = const., 
where o(2) = +(a- r(x, €) denotes 


the distance between the points x and ¢ in a plane 
section, and R(é) denotes the distance between 
the point € and the point of zero potential, the 
integration being taken over the complete tri- 
angular contour. Taking the locus of zero poten- 
tial at the center of one of the faces and setting 
the constant equal to — GQa?/6, the equation can 
be written in the form 


GO 

ai+a? 


Since the angle opening at the corners is 52/3 
the induced charge o(a2) must have singularities 


= 2/5 atz = aand x = 2a. 


Also, by symmetry, o(a) = o(2a — x). Because 
of the necessity of a considerable amount of ap- 
proximate integration in the present problem, it 
is convenient to define the singular approximation 
in such a way that it will vanish over a large part 
of the interval (0, 2a) but will have sufficient 
regularity in the open interval to enable the 


difference between the functions o(2) and s(x) to 
be approximated by a polynomial y(x). Such a 
function, having a continuous first derivative in 
the open interval, is obtained if we set 


8 7 
(=) + (4): () (4) 0 a/4, 
a 5 a § 
s(x) ={0, a/4Ka<7a/4, 
, 8 9 
(2—x/a)~*—-—(4) 
aj 


This function is formed, over the interval (0, a/4), 
by subtracting from the function (2x/a)~? the 
ordinates of the tangent line at x = a/4, and, 
over the interval (7a/4, 2a), by subtracting from 
the function (2 — x/a)-* the ordinates of the 
tangent line at x = 7a/4. 

Now, employing the notation of Problem 1, we 
set 


(3.20) y(x) = KiQ)yi, 
where = + K_,(A), = 2(2/a — 1), 
and ye, yi, and yo are the undetermined ordinates 
of y(x) at a = 0, a/2, and a, respectively, and 


require that c(7) = s(x) + y(a) = o(x). Then 
it follows that 

= —2@ {s(S)+y(S)} lo 

9 R(E =) dé, 


where, because of the restriction on the total 
charge, 


S odds = S + =0. 


Next we require that 


21) het 
Uk, 
2=—2@ j}s(c)+y(s)! log dé, 
k=1,. 
corresponding to equations (2.9), a the 


points a; are given, as before, by x = 0, a/8, a4, 

a/2, and 3a/4, respectively. The an 
be carried out for the most part by approximate 
methods, the work being arranged systematically 
in terms of matrices,!® yielding the system of 
equations 


(3.22) 
1.06226y2 — .79070y; — 3.48779yo + 2.98691% = 
50000 GOa, 
.99668y2 + .65565y, — 3.67312y¥9 + 1.22791ly = 
38281 GOa, 


> = 
ae 
° ‘ 
afl 
0 0 
4: 
O order — —— 
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68148y2 4- 1.35204y1 — 3.38659y0 + .65509y = 


.28125 GOa, 
.21117y2 + 1.30639y; — 2.05592yo + .21539y = 
.12500 G@a, 
03484y2 + .41273y; — .571338y0 + .04688y = 
03125 GOa, 
together with the condition 
23) 
S = .31111y2 + 1.42222y; + .26667y0 
+ .40626y = 0. 


Using (3.23) to eliminate y from (3.22) and 
applying the least squares method of Problem 1, 
we obtain the three equations 


(3.24) 
.40320y2 + 3.24045y; + .92688y> = 
— .09868 GOa, 
3.24045y2 + 39.29219y; + 24.8095lyo = 
— 2.22493 GOa, 
92688y2 + 24. + 33.56787y0 = 
— 2.63692 GOa, 


which give yz = .11808 G@a, y; = — .02757 GOa, 
yo = — .06144 GOa, and x = .04642 GOa, the 
value of x being obtained from (3.23). 

These values, in connection with (3.19) and 
(3.20), serve to determine the function (zx). 
Now, by approximate integration, it is an easy 
matter to determine the value of the function 9 
and hence of the stress function ® = 9 — GOr?/2 
at any point inside the contour. For since 9 = 
GQa?/6 at the center of a face, the value of 9 at 
any _ P is given by 

rp(§) 


—2¢3 


where rp(§) is the distance between the point P 
and the point &, and R(§) is the distance between 
the point € and the reference point at the center 
of a face. 

For example, the value of 9 at the center of the 
bar is given approximately by the integral 


9(C) GOa?/6 3(£) log V(E— a)? +a?/3 dt 


(§—a)?+a?/3 
This is evaluated easily by approximate methods 
with the help of matrices tabulated earlier in the 
problem’? and gives 9(C) = .22254 GO@a?. Since 
r=0 at this point, referring to (3.16), this 


gives directly the value of the stress function 
at the point. The exact value!’ is : GOna?, so that 


the error is less than 0.15%. 
Finally, we compute the value of the shear 
stress T, 


at the center of a face, where the maximum 
stress occurs. Now the integral 


P 


where P is a point on the boundary, §p(&) is the 
angle included between the normal at P and the 
line joining P and &, and Rp(&) is the distance 
between the points P and &, gives the mean value 
of the interior and exterior limits of the normal 
gradient of g at P.” Since the total jump at the 
boundary is given by 47x times the charge density 
at P,” it follows that the required interior limit 
is given by E(P) + 2x0(P), so that the shear 
stress t at P is given by 


<(P) = GOr(P) — E(P) — 2x(P). 
In particular, if Po is at the center of a face we 


have, because of symmetry and since 3 = 0 on 
that face, 


E(Py) = —4f 


bad 


sd5 


which gives, by approximate eee E(Po) = 
09914 G@a. Also, r(Po) = GOa/+/3 and (Po) = 
yo = — .06144 GOa, so that we obtain 7(Po) = 
£86424 GOa. The true value" is W3GQa/2 = 
.86603 GQa, the error thus being less than 0.21%. 

The exact solution to problem 2 can be found 
very easily by semi-inverse methods." However, 
the application of such methods is restricted to 
very special types of cross-sections. The present 
method obviously requires considerably more 
computation but may be of practical interest in 
cases where experimental results are not suffici- 
ently accurate, since it is a direct method and 
can be applied to any prismatical bar. 

4. Conclusion. The purpose of this paper is to 
present a practical method of obtaining numerical 
solutions to singular integral equations. We have 
not been concerned with existence and uniqueness 
theory, which has been well developed elsewhere. 

Although a certain amount of calculation is 
necessarily involved, the work can be organized 
and reduced principally to matrix multiplication 


\ 
re 
: d® dg 
= G@Or 
“ 
: 

‘ 
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and similar operations for which the modern 
computing machine is well adapted. Much 
simplification results from the use of the La- 
grangean form of the polynomial in connection 
with published tables of Lagrangean coefficients.! 

We have deliberately chosen problems for 
which the exact solutions are known. However, 
a large variety of important problems arising in 
static field theory, air-foil theory, elasticity, etc., 
which can be set up in a natural way in terms of 
singular integral equations, are at present 
unsolved. 

The principal practical difficulty remaining 
concerns the determination of the existence and 
nature of the singularities in the unknown 
function. In many cases such information can be 
obtained from physical considerations, as in 
Problems 1 and 2. Otherwise, one can assume 
the existence of a singularity of order ¢ at a 
certain point and determine ¢ by requiring that 
the integral equation be satisfied at that point, 
or that both sides of the equation have the same 
asymptotic behavior in the neighborhood of the 
point. 

No expression has as yet been devised for 
estimating the error inherent in the present 
method of approximation. Since the determina- 
tion of such an expression seems to involve the 
solution of a second integral equation, and since 
error estimates are generally of limited applica- 
bility, particularly when the solution is un- 
bounded or the data are empirical, a further 
study of the problem has not been made. How- 
ever, substitution of the computed values of the 
parameters in the m equations involved in the 
least square procedure gives an indication as to 
how well the integral equation is satisfied at the 
corresponding m points in the interval of inte- 
gration. If these points are chosen judiciously, a 
physical interpretation of the variation in the 
differences between the two sides of the equation 
will frequently indicate the degree of approxima- 
tion obtained. 

Finally, we point out the fact that many 
singular integral equations in which the range of 
integration is infinite can be transformed to 
types to which the methods outlined will apply. 
For example, to determine the induced charge 
distribution on an infinite half-plate it is necessary 
to solve the integral equation 


o(x) = —2 S o(E) log 


(where the locus of zero potential is taken to be 
the line x= 1). If we set § =7/(2Qa—7), x= 


t/(2a — t), and «(t) = (2a — , the 


equation takes the form 


t 2a 
= —4a log 
«(z) log 


4a log (2—t/a) K dz, 


dz + 


the approximate solution to which can be readily 
obtained with the use of the data computed for 
Problem 1. 
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